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ABSTRACT 


This investigation presents a Spectral Method for solving 
the one-dimensional fallow-water wave equations. The spectral 
method is based on the Chebyshev Collocation technique and a 
first-order finite difference time stepping. The method developed in 
this study is applied to route a Log-Pearson type III hydrograph 
through a wide rectangular channel. The spectral solutions are 
compared with the solutions obtained using the Preissmann finite 
difference scheme which is first-order accurate in space and 
second-order accurate in time. An error analysis is made through a 
parametric study. It is concluded from this study that the spectral 
method performs better than the Preissirann scheme as long as the 
time stepping errors are kept low. However, a second-order 


t ime -accurate 

finite difference scheme which may have 

only 

a 

first-order accuracy in space 

, performs better 

than 

a 

spect 

ral 

method which 

uses first-order 

finite difference 

t ime 

stepping 

i n 

cases of fast 

rising f 1 oods and 

friction dominated 

flows 






IV 


ACKNOWLEDGEMENT 

I express my heartfelt thanks and gratitude to my guides 
Dr. B.S. Murty and Dr. V. Eswaran for their guidance, support and 
encouragement during my thesis work. They have always been ready to 
help me. I cannot forget the patient and sympathetic hearing that I 
got from them whenever I had any problems - personal or related to 
work. I am also grateful to all faculty of I IT Kanpur. 

I thank my friends with whom I shared my enjoyments and 
worries. In particular, Bimal Modi, Ashok Keshri and V. Jagid have 
been a great help in preparation of this thesis. 

Lastly, I convey my love and regards to my parents for 
their constant inspiration and encouragement. 


Jiweshwar Sinha 



D^dicat€^d 


to th& 


UNEMPLOYMENT 



TABLE OF CONTENTS 


Page 

CERTIFICATE i i 

ABSTRACT i i i 

ACKNOWLEDGEMENT i v 

LIST OF FIGURES v i i i 

LIST OF SYMBOLS ix 

CHAPTER I INTRODUCTION 1 

CHAPTER II GOVERNING EQUATIONS AND THE FINITE DIFFERENCE 

2.1 St. Venant Equations 8 

2.2 Non-Dimens ional izat ion 10 

2.5 Preissmann Scheme 12 

2.5.1 Boundary Conditions 15 

2.5.2 Algorithm 16 

CHAPTER III SPECTRAL METHODS 

5.1 The Collocation Method 20 

5.2 Chebyshev Polynomial and Evaluation 

of Derivatives 21 

5.5 The Chebyshev-Col locat ion Method with 

F.D. Time Stepping 24 

5.4 Preconditioned Residual Minimisation 26 

CHAPTER IV SPECTRAL SOLUTION OF THE GOVERNING EQUATIONS 

4.1 Residual Determination 28 

4.2 Finite Difference Preconditioning 50 

4.5 Preconditioned Residual Minimization 55 



VI 1 


CHAPTER V RESULTS AND DISCUSSION 

5.1 Non-Dimensional izat ion of Equations 56 

5.2 Weighting Parameters a, ft and 9 58 

5.5 Results from Equations without 

Source Term 59 

5.4 Results from Equations with 

Source Term 47 

5.5 Computational Time 55 

5.6 Conclusions 55 

CHAPTER VI SUMMARY 56 

REFERENCES 58 

APPENDIX A 66 


APPENDIX B 


68 



LIST OF FIGURES 


Figure Page 

Z.l Definition sketch 9 

2.2 Computational grid 14 

4.1 Flow chart of the algorithm 54 

5.1 "Exact solution" for hydrograph-l 41 

5.2 "Exact solution” for hydrograph-2 42 

5.5 Spatial resolution error in depth 

computation, hydrograph-l 44 

5.4 Spatial resolution error in velocity 

computation, hydrograph-l 45 

5.5 Spatial resolution error in depth 

computation, hydrograph-2 46 

5.6 Time-stepping error in depth 

computation, hydrograph-l 48 

5.7 Time-stepping error in depth 

computation, hydrograph-2 49 

5.8 "Exact solution" for hydrograph-l 

source term included 51 

5.9 Spatial resolution error for hydrograph-l 

source term included 52 

5.10 Time-stepping error for hydrograph-l 

source term included 54 


LIST OF TABLES 

4.1 Coefficients of possible sets of 


non-dimens i onai ized St. Venant equations 


57 



LIST OF SYMBOLS 


General Symbols 

h flow depth 

V flow velocity 

t computational time 

X co-ordinate in the flow direction 

n Manning's coefficient of roughness 

L channel length 

Ax non-uniform grid separation 

At time increment 

N number of collocation or grid points 

N number of spectral coefficients or frequency modes 

c 

J Jacobian 

Af correction to solution 

R residual 

Chebyshev polynomial of the first kind 
c, ,c specified constants 

k m 

r known quantities on R.H.S. 

a-d known quantities on L.H.5. 

Special Symbols 

a -g non-dimensional constants 

o o 

p dimensional constant 

o 

S non-dimens ionl ized source 

r 

t ratio of the time at C.G. 

r 

at peak discharge 

t^ time to peak discharge of 

q peak discharge 

p 


term 

of hydrograph to the time 

inflow hydrograph 



Greek Symbol 


a. 

ft 

0 

4 > 

b) 


weighting factor 
weighting factor 
weighting factor 
angle of channel 
minimum residual 


in F.D. 

in F.D. preconditionar 
in spectral method 
-axis with horizontal 
relaxation factor 


Subscripted Symbol 
sp spectral 

j grid point 

p time level 

t partial derivative with respect 

X partial derivative with respect 

o initial uniform condition 

e exact solution 


to t 
to X 


Supercript Symbol 

n number of iteration 

u upstream equation 

d downstream equation 

c continuity equation 

m momentum equation 

non-dimensional form 


first derivative 



CHAPTER I 


INTR<»UCTION 

Unsteady free-surface flow has been the subject of 
extensive studies for many decades because of its practical 
importance. Unsteady flow may be produced in a channel by intense 
rainfall (floods), by the failure of a control structures 
(dambreak wave), by the movement of control structures <Ex: sluice 
qates) in irrigation canal networks or by other natural or 
man-made disasters. The analysis of unsteady flow becomes 
important for determining the limits of potential damage due to 
floods, for proper design and operation of various hydraulic 
structures and for determining the morphological changes in the 
rivers. An unsteady flow analysis of a flood basically involves 
the determination of hydrograph at any cross-section in the 
channel if the hydrograph is given at an upstream location. Over 
the years, this analysis has progressed from an idealized 
analytical solution of simplified equations to the numerical 
models for complex situations. The numerical models involve the 
solution of partial differential equations describing the unsteady 
flow for prescribed boundary conditions. Depending on the 
idealization of the physical situation under study, one or more 
terms in the complete governing equations may be neglected. 

The basic theory of unsteady flow in open-channels was 
first formulated by Barr’e de Saint Venant in 1871. He derived the 
governing equations for one-dimensional unsteady flow in prismatic 
channels. These equations, also known as shallow water wave 
equations, describe the conservation of mass and momentum and are 
based on the assumption of hydrostatic pressure distribution in 



the vertical direction. The two governing equations are non-linear 
and are hyperbolic in nature. The earliest attempt at finding 
solutions to these equations was made by Massau in 1889 <Mahmood 
and Yevjevich, 1975) using graphical integration of the 

characteristic equations. Ritter <1892) used these equations to 
analytically determine the dam-break wave movement in a 
frictionless. horizontal, prismatic channel. Later, Dressier 
<1952) extended the Ritter solution to include the friction 

effect. However, the analytical solutions are available only for 

highly idealized cases and simple boundary conditions. The first 

engineering application of shallow water equations had to wait for 
the development of digital computers so that numerical techniques 
could be applied for solving them. The first attempt to 
numerically solve the one-dimensional St. Venant equations for 
practical problems was made by Stoker <1955) and Issacson, et al . 
<1954). They refined the finite difference approach of Ihonrtas 
<1957) and applied it for routing floods in Ohio and Mississippi 
rivers. Since then, much effort has been expended in the 
development of the numerical techniques for solving the shallow 
water wave equations. These methods are well described in many 
text books <Mehmood and Yevjevich. 1975, Abbott, 1979, Cunge , et 
al . , 1980 and Chaudhry , 1987) and only a brief review of 

literature is given below. 

The characteristic method of solution was a natural 
choice for the early researchers because of the hyperbolic nature 
of the equations. One-dimensional explicit characteristic models 
were developed by Liggett and Moolhiser <1967), Streeter and Mylie 
<1967) and Ellis <1970) while the implicit characteristic models 
were developed by Amein <1966) and Wylie <1970). The 


characteristic methods require either a curvilinear qrid or 
interpolations in the x-t solution domain because of varying wave 
speed. Therefore, their application for flood routing is not 
popular in recent times. The characteristic method also requires 
special treatment in case shocks such as bores are present in the 
solution. 

The shallow water wave equations are similar to the 
equations of gas dynamics. Therefore. the finite difference 
techniques developed in gas dynamics can be used for solving the 
flood routing problem also. This was first recognized by Stoker 
and Issacson, et al . . In a finite difference scheme, the governing 
equations are converted into algebraic equations by replacing the 
part ial - di f f erent ial terms by their finite difference analogs. 
These algebraic equations are then used to advance the solution 
from a known time-level to an unknown time-level. The finite 
difference scheme can be either explicit or implicit in nature. In 
an explicit finite difference scheme, the spatial derivatives and 
the coefficients are evaluated at the known time-level. Tollowing 
the pioneering work of Lax <1954) and Lax and Wendroff <1960), 
several explicit finite difference models were developed for 
shallow water wave equations. Notable among them are those of 
Liggett and Moolhiser <1967), Martin and DeFazio <1969), Dronkers 
<1969). Strelkoff <1970) and Jhonson <1974). Many of these schemes 
use either the first-order time accurate Lax diffusive scheme or 
the second-order time accurate Leap-frog scheme <Liggett and 
Cunge, 1975). Recently, Fennema and Chaudhry <1986) introduced 
MacCormack <1969) and Gabutti <1985) schemes for the analysis of 
unsteady free-surface flows with shocks. Both MacCormack and 
Gabutti schemes are second-order accurate in time and space. All 
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the above explicit models have a restriction on the size of the 
computational time step because of numerical stability. Therefore, 
they are not suitable for application to routing flood waves, 
having a duration in order of days. This necessitated the 
development of implicit models. 

In an implicit scheme, the solution is advanced from one 
time-line to the next simultaneously for all points along the 
time-line due to the fact that the spatial derivatives and the 
coefficients are evaluated at the unknown time-level. The implicit 
schemes are generally unconditionally stable and therefore, larger 
computational time steps can be used. The use of implicit models 
was first suggested by Issacson, et al. <1956) and became popular 
after the development of Preissmann four-point implicit scheme 
<1961). Later, Vasil iev, et al . (1965) and Abbott and lonescu 
(1967) developed six-point implicit schemes. However, these 
schemes could not become popular because of the requirement of 
uniform finite difference grid. Also, the inclusion of boundary 
conditions is more complicated in the six-point schemes. Amein and 
Fang (1970) are the first to use the Newton-Raphson technique to 
solve the system of non-linear equations in an implicit model. 
Recently. Fennema and Chaudhry (1990) used the Beam and Warming 
implicit scheme (1976) for solving the unsteady free-surface flow 
problems. Notable among the other implicit models are those 
developed by Baltzer and Lai (1968), Tread (1975.1976) and Chen 
and Simons (1975). 

The numerical stability and accuracy of various implicit 
schemes have been studied by Cunge (1966), Abbott and loneacu 
(1967). Tread (1974) and Liggett. et al . (1975). Although, 
implicit schemes are usually unconditionally stable, Abbott and 



lonescu <1967) have demonstrated the danger of introduction of 
significant numerical errors when excessively large time steps are 
taken. A means of circumventing these errors is also described by 
them. Chaudhry and Contractor {1975> and Fread (1974) have shown 
that large computational time steps can also result in numerical 
instability. This is especially so when the spatial derivative 
terms are not sufficiently weighted towards the future time-line 
and the changes are rapid. Fread <1974) investigated the 
convergence properties of various four-point implicit schemes by 
determining the truncation error. The backward implicit scheme has 
a first-order truncation error where as the box implicit scheme 
has a second-order truncation error. 

Taylor, et al . <1974) were probably the first to attempt 
at solving the shallow water equations using finite element 
technique. Later, Cooley and Moin <1976), King <1976), Keuning 
<1976) and Meissner <1978) used the finite element method based on 
Galerkin approach. It was concluded from their studies that the 
results obtained by the Finite difference methods are, in general, 
superior to those obtained using the finite element models. 
Katopodes < 1 978 , 1 980 , 1 984a , b) addressed the basic difficulties 
encountered in the application of finite element methods for 
solving the shallow-water wave equations and proposed a new 
dissipative Galerkin scheme for open- channel flow. However, the 
finite element models are yet to become as popular as finite 
difference models at least for one-dimensional computations. They, 
however, are preferred by The U.S. Army Corps of Engineers <1987) 
in two-dimensional modeling with complex boundaries. 

Both finite difference and finite element methods have 
only a finite-order accuracy. A considerable increase in both 


algorithmic and progratttmi ng complexity is involved if the order of 
accuracy of these methods is increased. This, however is not the 
case with the Spectral Method. Claison, et al . <1970) and Gottlieb 
and Orszag (1977) laid the foundation of the spectral method and 
developed the transform method and the spectral collocation 
techniques. The book on spectral methods by Canute, et al . <1987) 
traces the history of the methods and discusses in detail the 
various spectral techniques starting from first principles. The 
book also presents several application of spectral methods to 
fluid flow problems. Only a very brief review relevant to the 
present work is given below. 

Orszag <1980) outlines the spectral collocation 
technique using Chebyshev polynomials and applied it to 
non-linear, non-constant coefficient problems. Spectral accuracy 
is achieved only in space while the time differencing is done by 
using the standard implicit finite difference schemes. Orszag 
<1980) also introduced the concept of preconditioning For 
efficiently solving the full matrix equations resulting from 
spectral approximation. It was shown that the extra work required 
in solving the spectral equations above that of solving the lowest 
order finite difference approximation is only 0<N log N) where N 
is the number of grid points m the numerical scheme. Moin and Kim 
<1980) used a pseudo- spectral formulation based on Fourier series 
and Chebyshev polynomial expansions and solved some channel flow 
problems. A mention was made about the aliasing error but it was 
suggested that this error is limited due to the damping of the 
physical system. 

In the studies reported herein, an attempt is made to 
develop a Spectral Method for the solution 


of one- di mens ional 



shal low- water wave equations. This spectral method is based on 
Chebyshev Collocation technique and finite difference time 
stepping. The spectral solutions for a hypothetical flood routing 
problem are compared with the finite difference solutions obtained 
using the Preissmann scheme. The spatial resolution and the time 
stepping errors in numerical solutions are analysed through a 
parametric study. In this report, the governing shallow-water wave 
equations and the Preissmann scheme are presented in Chapter II. 
In Chapters III and IV, details of Spectral Method using Chebyshev 
Collocation technique are presented. In Chapter V, the numerical 
results are discussed and an error analysis is made. 



CHAPTER II 


GOVtIRNING EQUATIONS AND THE FINITE DIFFERENCE 


2*1 SI. Venani Equations 

The gradually varied unsteady Flow in an open channel is 
described by the Saint Venant equations. Considering the unsteady 
flow in a wide rectangular prismatic channel as shown in fig 2.1, 
these equations are derived using the following assumptions: 

<a> the flow is one-dimensional in the direction of 

the channel; that is, the velocity is uniform over a 
cross* sect i on and the water level across a section is 
hor i zontal , 

<b> vertical accelerations are negligible, 

(c) steady state resistance laws may be used to model 

unsteady effects, 

(d> the average bed slope, of the channel is small such 
that the cosine of the angle made by the averaged bed 
profile with the horizontal can be regarded as unity. 

By making the above assumptions, the Saint V/enant 
equations are obtained by applying the principles of conservation 
of mass and conservation of momentum to a differential control 
volume (Chaudhry, 1986). The resulting equations in a 
non- conservat ion form are: 


dh 

at 


+ v 


dh 

dx 


+ h 


dv 

dx 
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(2. 1 .la> 


and 


dv 
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+ Of 


dh 

dx 


= q <S, 




dt 


+ V 


dx 


(2.1 .lb) 
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where h is the depth of flow; v is the flow velocity; q is the 
acceleration due to gravity; is the slope of the 
channel- bottom: is the slope of the energy grade line; x is the 
distance measured positive in the downstream direction, and t is 
the time of computation. In the present work, is evaluated 
using Manning's formula for a wide rectangular channel. 



where n is the Manning's roughness coefficient. 


< 2 . 1 . 2 ) 


2 Non-Di mens! onal 1 :zat.l on 

The Saint Venant equations are non-dimensional ized in 

the present study using channel length, as distance scale 

parameter and the initial steady uniform flow depth, as the 

depth scale parameter. The initial steady state uniform flow 

velocity, V is taken as the velocity scale parameter. V is 

o o 

obtained using the Manning's equation, 

. 2/3 i/2 

V = -i_ H S <2.2. la) 

o n o o 




t 


IZ.Z.Zh) 


= t / T, 


h - h / 


<2. 2. 2c) 



<2.2.2d> 


Substitution of equations <2.2.1) and <2.2.2) in equations <2.1.1) 
and <2.1.2) result in the followinq general non-dimensional St. 
Venant equations 


and 


dh , . “ dh , r dv 

+ b V “t c h 

9t ® dx * dx 


<2 .2. 5a) 


dv . - dv 

— _ + e V — ^ 

Ot ** dx 


+ f 


ah 

ax 


= q <S 

o o 


o f 


<2*2. 5b> 


where a . b • c , d , e , F and q are dimensionless constants, 
o o o o o o o 

but is dimensional constant and it becomes dimensionless when 

combined with S, where 
t 


n 


2 -2! 
V 


<2.2.4) 


Depending on the v\iay the non- dimensional ization is 

carried out, as many as 12 variations of equation (2.2.5) are 

possible. The best set of equations among these is considered in 

the present study. The basis of selection of this set is discussed 

in Chapter V. Only the definitions of coefficients a - q and p_ 

o o o 

corresponding to the best form of non-dimensional St. Venant 
equations are given beloM. 
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= 1.0 


b = V T / L 

O O O ' o 


c = b 

o o 


d = 


eo 


V / qT 

Vo / qLo 

H / L 
o' o 


q =1.0 


C2.2.^a> 


<2.2.5b> 


<2.2.^c> 


<2.2. 5d> 


<2.2.5e> 


<2 .2. 5F> 


<2.2. >q> 


and 


p = V* / (2.2. >h) 

o o' o 

The riqht hand term of the equation <2.2.5b> is a source term 
which is highly non-linear in nature and requires a different 
treatment. Therefore, it is represented by a special symbol . 

S = q <S - p S, ) (2.2.6) 

r -O O O f 

2* 3 HreJ ssmann Sfcheme 

The governing equations presented in the previous 
section constitute a set of non-linear, first-order, hyperbolic 
partial differential equations, analytical solutions for which are 
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available only for highly idealized cases. ThereFore. they are 
generally solved using numerical techniques. In this section, one 
such scheme, originally developed by Preissmann <1961). is 
presented . 

Preissmann scheme (1961) is basically a four-point 
weighted implicit finite difference scheme. It is used to 
transform the non-linear partial differential equations of St. 
Venant into non-linear algebraic equations. Fleferring to fig. 2.2, 
a function f(x.t) is represented by its values on a arid with F^ 

j 

denoting the value at the jth grid point and the pth time level. 
Ax. IS the non-dimensional grid spacing between the <i41>th and 
the ith point. Then for a non-dimensional time step At, the 
Four- point weighted difFerence approximations For the time and 
space derivatives f^ and f for a "cell” defined by the Four 
points <p,j>, <pHl,j>, <p,}fl>. (p^l,j^l> are 


F 


t 


f - f ** + f - f ** 
» j j 

2 At 


<2.5. la> 


f 

X 


oc 




<i-«> 



<2.5. lb> 


and 


fP"** 4- fP'** ^P + f P 

f(x.t> = « J — i — -J < <l-«) ^ ^ <2.5.1c> 

where oc is a weighting factor varying from 0.5 to 1. 

A scheme using an oc value of 0.5 is known as the ’box 
scheme' while oc - 1 results in the ’fully implicit scheme’, for oc 
- 0.5, the scheme is of second-order accurate in time, otherwise 





it remains first-order accurate. Tread <1974> recommends <x values 


nearer to 0.^ to insure unconditional linear numerical stability 
and also provide good accuracy. 


2. 3. 1 Boundary Conditions 

The four-point implicit scheme presented in the previous 
section is based on forward finite differencing. Therefore, from N 
nodes only (2N-2) algebraic eguations are obtained using the St. 
Venant equations <2.2. 5>. However, the number of unknowns at the 
advanced time level are 2N. The additional two equations required 
for completely solving the system are obtained from boundary 
conditions. The specification of boundary conditions depends on 
whether the flow is subcritical or supercritical. In the present 
work, we are interesed only in the case of subcritical flows. 
Therefore, one boundary condition is specified at the upstream end 
and another is specified at the downstraem end <Stoker 1957). The 
upstream boundary condition is specified by a hydrograph 
(typically of type III) and a second order extrapolation d h / dx' 
- 0 IS applied at the downstream end. 

The general form of the type III hydrograph is 
represented by 


q(t) -- q^ 


l-i<q -l)(t/t 
P P 


,Cl/<t -1)3 

) r 
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Cl 


<t/t )3C1 
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■<t 




< 2 . 5 . 1 . 1 ) 


in which q<t) is the discharge per unit width at time t; q^ is the 
initial discharge; t is the time at which peak discharge occurs; 

p 

q is the ratio of discharge at peak to initial discharge and t 

r T 

is the ratio of time at Centre of Gravity of the hydrograph to the 





time of peak discharqe. 

The second order extrapolation / dx* - 0) at the 

downstream end results in the following equation. 


(X - 
N 


N- 


)h 

1 N- 2 


<X 
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X )h 
N-2 M- 


<x - X >h =• 

N-l N-2 N 


<? .5. 1 .2) 


where subcript represents the grid-point. 

It should be noted here that the above second- order 
extrapolation at the downstream end is strictly not valid. 
However, it gives proper results as long as we are not interested 
in the upstream effects of the downstream boundary. Therefore, the 
above methodology can be used only for routing floods in channel 
reaches which are away from dams, weirs, channel junctions etc. In 
such cases, a rating curve at the downstream end can be used as a 
boundary condition. 


2*3*2 Algoritrhm 

Equations <2.2.5> are discretized using the equations 
<2.?.1>. This results in a system (2N-2) of non-linear algebraic 
equations corresponding to N grid points. Two separate non-linear 
algebraic equations are obtained from upstream and downstream 
boundary conditions (2. 5. 1.1) and <2.5.I.2>. The system of 
equations can be represented in a functional form as given below, 

B <h , V ) = 0 C2.5.2.1a> 

1 11 

= 0 . 1 < i < N-l <2.3.2.1b> 

= 0 , 


C. (h.,v.,h. ,v. > 

} 3 3 j+i J+4 


M. <h. , V. , h. , V. ) 
j 3 3 j + i 


1 < j < N-l 


<2.5.2.1c> 
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and 




N- 2 


N- 1 


h > = 

N 


<2 . 3.2 . Id) 


where subscript denotes the grid-point and the barred notation is 
for the normalized value of the variable. In <2. 3. 2.1) C and M 
correspond to the functional form of the discretized St. Venant 
equations and B corresponds to the functional form of boundary 
conditions. Complete details of these equations are presented in 
the Appendix A. 

The non-linear equations <2. 3. 2.1) are solved in this 
study using Newton-Raphson method. The Newton- Raphson method is 
basically a functional iterative technique to solve a system of 
non-linear equations. This technique is derived from a Taylor 
series expansion of the non-linear function in which all terms of 
second and higher order are neglected. The resulting algorithm is. 


3 <f’^) Af = - F<f”> 


< 2 . 5 . 2 . 2 ) 


in which f’^ is a vector quantity, 3 is the Jacobian <a coefficient 
matrix made up of the partial derivatives evaluated with f"^ 
values). F<f'^> is the non-linear vector equation evaluated with f” 
values, and Af is a vector containing the 2N unknowns. Equation 
<2. 3. 2.2) represents a system of equations in which the unknown 
vector Af is linear . Thus, the Newton-Raphson method generates a 
2N X 2N matrix equation. This matrix is banded with most of the 
elements being zero except for four elements in each row along the 
main diagonal of the matrix. In the present work, application of 



the Newton-Raphson method to the equation <2. 5.2.1) results in the 
follovsfing matrix equation. 
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^2 .5. 2 .5) 


Tn the above equation, coefficients a'^ . , c* and d\' are partial 

j j j J 

differentials of the discretired continuity equation (2. 5. 2. lb) 
with respect to hf ^ , h*^'* * and * respectively; 

1 i j ■' 1 j ' 1 

coefficients a"*, b”*, c"* and d”* are partial differentials of 

j j j j 

normalized dynamic equation <2.3.2. Ic) with respect to h*^^ * , v^^ ^ , 

J J 

hf*^ and v*?** respectively; coefficients a^ and are partial 

j+i j+i 1 1 

differentials of the normalized upstream boundary condition 
<2. 3-2. la) with respect to and respectively; 

dt cl cl 

coefficients a , a and a are partial differentials of 

N-2 N-1 N 

the normalized downstream boundary condition <2.3.2. Id) with 

i j. 1 P+l i_P+* j i-P+i i • I V c m _i d 

respect to h*^ , h and h*^ respectively, r , r., r. and r are 

the residuals obtained from upstream boundary condition 

<2. 3. 2. la), continuity equation <2. 3. 2. lb), momentum equation 
<2. 3. 2.1c) and downtream boundary condition <2. 3. 2. Id) at time 
level <p+l> respectively. Complete details of the coefficients are 


presented in the Appendix B. 





The solution of St. Venant equations by a finite 
difference scheme essentially involves determining the depth and 
velocity values at all the grid points at the unknown, (p+1) time 
level using the depth and velocity values at the grid points at 
the known, <p> time level. This is accomplished using the 


following procedure. 


<a> 

<b) 


Assume values of h^^* . for all j = 1 to N. 

j j 


Use 

h*^, 

J 

v*^, h*^"" and v*^"" to 

4 J J 

determine r'^ , r® , 

4 J 

m 

r . , 
a 

d 

and 

the 

elements of the 

coefficient matrix 

in 

the 


<c) 


(d> 


<e) 


<f > 


equation <2.5.Z.5>. 

Determine Ah, and Av . for all j = 1 to N by solving the 

J J 

matrix equation <2.3.2.3>. A standard library subroutine 
is used for this purpose. 

The corrections Ah. and Av . are added to hf** and v*^*^ 

J J J J 

to get the new guess values for and v*?"*"*. This 

J J 

entire process is repeated till convergence. 

The values of and obtained after the 

J J 

convergence are referred to as the finite difference 
solutions of the St. Venant equations. 

Once the solution is obtained at the (p+l> time level, 
it is advanced to the next time level and the procedure 
is repeated till the end of simulation period is 
reached. 


The finite difference solution at any time forms the 
initial guess for the Spectral solution at that time. The basics 
of Spectral method are presented in Chapter III and the complete 
algorithm is discussed in Chapter IV. 



CHAPTER III 


SPECTRAL METHOD 

This chapter contains the details of the spectral method 
including Chebyshev polynomials. the Chebyshev-Col locat ion 
techniques and the Preconditioned Residual Minimization method. 

3. 1 The Collocation Method 

The spectral collocation technique is generally used for 
solving fluid flow problems. The technique uses expansion 

functions such as orthogonal Chebyshev and Legendre' polynomials, 
and treats the function values at certain physical points as the 
fundamental representation of the solution. 

In the Collocation approach, the truncated series of the 
solution is represented as; 

N -1 

h <x) = r h, B, <x> , -1 < X < 1 <5.1.1) 

N , “ k k 

k=0 

where h is the spectral representation of h, h are the spectral 

IhT jC 

coefficients, B. are the orthogonal expansion functions and N is 

k C 

the number of terms in the expansion. 

For any differential equation 

M<h) = f<x> , -1 < X < 1 <5.1.2) 

where M<h> is a differential operator which includes spatial 
derivatives of h, and f is the known forcing function. If the 
approximate solution <5.1.1) does not satisfy <5.1.2) everywhere. 


then the residual 
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R <x> = f (x) - M <h ) (j 1.5> 

does not vanish everywhere. In (5.1.5> f and M are the spectral 
approximations of f(x> and M(h) respectively. For the collocation 
method, the prime requirement is that the P.D.E. be satisfied 
exactly at each of the collocation points , i.e. 

* = N ‘’-I-*’ 

J 

The spectral coefficients hj^ can be found either by 
analytic integration using the orthogonal property of the basis 
functions or by finite transformation which relate the 
coefficients to the function values at collocation points. The 
first approach is called the ” Galerkin ” and the second the " 
Collocation " . 


3.2 CKdbyshev Polynomials and Evaluat-lon of Derivatives 

The Chebyshev polynomials are a family of orthogonal 
polynomials in the interval <-l,l). Chebyshev polynomials satisfy 
the orthogonal relationship. 


1 

f 

-i. 


T, <x) T (x) 
k m 


dx = 


n 


/ 


C 6 
k km 


where C, = 
k 


1 - x" 




2 if k = 0 


(3.2.1) 


otherwi se 


where T (x> are the Chebyshev polynomials of the first kind. If 
k 

the approximate solution h (x) is represented as 


h <x> 

N 


i K 

k=0 


(5.2.2) 
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then using (3.2.1) the equation for the Chebyshev coefficients 
becomes 




h 


1 

c 

m 


1 h (X) T <x) 

/ —ii 12 

y 1 - X 


(3.2.3) 


The Gauss-Lobatto quadrature rule which provides a means 
of representing the continuous integral of a function in terms of 
a discrete summation, gives the Chebyshev-Gauss-Lobatto 
Collocation points 


X . = cos 
j 


[-V-] 


(3.2.4) 


Thus the spectral coefficients are determined by the forward 
transformat i on ; 


N 


hu = S . h . 


k = 0,1,2 N -1 

c 


(3.2.5) 


where h. are the values of the function at the collocation points 

j 

X . and 

J 


k j 


(N-1) C. C, 

4 k 


r Tr( j-l>k T 

'=°H N-i" J 


j = 1,2 N 


(3.2.6) 


k = 0,1,2 N -1 

c 


where 


p2j = l,N = 

= i * C, = ^ 

■’ L 1 2 < k < N-1 1-1 2 < k < 


2 k = 1,N -1 
c 


(3.2.7) 


N -2 
c 


The grid point values of the function can be obtained 
from the spectral coefficients by the inverse transformation 


given by 



h. 

J 


N - 1 
c 


1 2 

*«*-»•* *1 


(3.2.8) 


= E 

k=0 






i 


N 


where the inverse transformation matrix is defined as 



r n<j-l)k 1 

N-1 I 


j = 1.2 N 

k = 0,1.2 N^-1 


(5.2.9) 


Both the forward and the inverse transforms can be 
implemented by using either Matrix Multiplication or the Fast 
Fourier Transform depending upon the number of coefficients and 
the number of grid points. 

The derivatives of the function h having Chebyshev 
polynomials can be written as 


^ T^(x) 

k=0 


(3.2.10) 


< 1 > 

where hj^ are the spectral coefficients of the derivative. 
From the relationship. 


2 T, <x> 

k 


1 

k+1 


T (x) 

k + l 


1 

in 


T 


k-4 


<X) 


k > 1 (3.2.11) 


we can obtain the recurrence relationship 


2 k h, = C 


'< ± > 
'k-l 


- h 


< ± > 
k4-i 


k > 1 


( 3 . 2 . 12 ) 


This relation gives an efficient way of differentiating a 

I ± y 

polynomial of degree N in Chebyshev space. Since h^^ =0 for 

k > N^-1 , the non-zero coefficient are computed in decreasing 


order by the recurrence relation; 
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C +2 <k + l) h. , 0<k<N-l (5.2.15) 

k k k+2 k+1 C 

After knowing the coefficients of h, the coefficients of 
the derivative of the function can be known using (5.2.15). These 
coefficients are then transformed back into physical space using 
an Inverse Transform to obtain the grid point values of the 
derivative. By this way a very accurate approx ittiat i on of the 
derivatives of the function is determined. 

3.3 The Chebyshev-Col location Method with F. D. Time Stepping 

Let us apply the Chebyshev-Col locat i on method to the 

equation 

h^ = M(h) , -1 < X < 1 (5.5.1) 

where M(h) is the differential operator containing spatial 
derivatives of h(x,t) and h^ is the partial derivatives of h with 
respect to time. 

The Chebyshev expansion of a function h(x,t) in the 
interval (-1,1) is given by 

00 ^ 

h(x,t) = r h, (t) T (x) (5.5.2) 

, k k 

k=0 

If the function h(x,t) is infinitely differentiable on (-1.1), 
then the Chebyshev coefficients decay faster than algebraically 
for increasing k. This ensures the spectral accuracy of the 
method . 

In practice, the truncated series of the solution is 
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represented as 

N - 1 

h„<x.t> = E h, <t) T, <x) . -1 < X < 1 <5.3.5) 

k=o ^ 

where N is the number of modes, 
c 

In solving <5.5.1) by a spectral method, the time 
derivative is discretized by finite difference technique. For a 
simple implicit time-stepping scheme, <3.5.1) can be 
approximated as 

h*"*^ = h*’ + At f M (h*’’^^) 1. , j = 1.2 N <5.3.4) 

j j L «p N Jj 

where the subscript j refers to the values at the collocation 

points x^ and the superscript p refers to the time-level. The 

number of points N is equal to the number of coefficients . The 
operation M <h ) is evaluated spectrally at the collocation 

op N 

points by using the collocation method explained in the previous 
sect i on . 

The scheme <5.5.4) requires an iterative process for 
obtaining the solution because of the unknown presents in the 

spectral operator M . However the spectral operator M is 

op op 

generally difficult and expansive to invert. This difficulty is 
circumvented by adopting finite difference approximation of the 
spectral operator for use in the iterative solution of <5.5.4). 
The approximate operator is chosen such that it is close to the 
original spectral operator and is inexpensive and easy to invert. 
This is known as Preconditioning. This preconditioning reduces the 
computational cost of the iterative solution procedure. 
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3.4 Preconditioned Residual Minimization 

Preconditioned residual minimization is the basis of the 
iterative algorithm used in the present work. This technique is 
briefly explained here. Let us consider a spectral equation as : 

M <h> = f <3.4.1) 

mp 

where M is the spectral representation of a differential 

« p 

operator M and f is a known forcing function. 

For any guess solution h” the residual in the equation 
is determined as 

R” = f - M <h") <3.4.2) 

The next step is to evaluate the corrections Ah’^ to be applied to 
the solution. This is done by using a finite difference solver. 
To ensure rapid minimization of the residual and to accelerate the 
iterative process, a relaxation factor co is incorporated in the 
solution as 

= h" + w" Ah’^ <3.4.5) 

In the collocation method, the minimum residual method 
is the process of the minimization of the sum of squares of the 
residuals at the collocation points, this sum is called the "inner 
product” and is shown as <R, R) . If R. and R. are the residuals 

j j 

at the collocation points obtained through the solutions of h*^ and 
h’^ + Ah’^ respectively, we choose such that <r 7** R*^**) is 

minimized. For a linear operator M 





< 5 . 4 . 4 ) 


= r"- ^ <AR) 


where AR is determined from 


AR = R - R’ 


<3. 4.5) 


The inner product of <3.4.3> can be given as 


[R?"^ R*^"" ) = (R, R*^) - 1 (AR, R ) + w* (AR, AR) <3.4.6) 


The magnitude of (3-4.6) is minimum for a value of the relaxation 
coefficient 


.n <R , AR) 

* <AR, AR) 


<3.4.7) 


Using and <3.4.3), the next guess h’^** is obtained 
and the process is repeated until the square root of the sum of 


the residuals decreases no further. 



CHAPTER IV 


THE SPECTRAL SOLUTIOM OF EQUATIONS 

This chapter presents the salient features of the 
algorithm of the Spectral Method used in the present Mork. The 
algorithm is based on Preconditioned Residual Minimization. The 
topics covered include the Residual determination. Finite 
Difference Preconditioning and Preconditioned Residual 
Minimization. At the end of the chapter, a flow chart of the 
algorithm used is included. For convenience, the bar has been 
dropped from the barred variables in Chapter II without changing 
the nature of the variables. 


4. 1 Residual Determination 


The governing equations are given as below. 


dh . . i. dv , ^h _ 

a "jET* f h h jk t c V Ts — — 0 

o dt o dx o dx 


. dv , dv , dh C - n 

d -37- + e V s- + f -s— - S =0 

O dt O dx O dx r 


<4.1 .la) 

(4.1 .lb> 


As mentioned earlier, the spectral scheme is applied by 
discretizing the temporal derivative with finite difference 
approximation and the spatial derivatives are estimated using the 
spectral Chebyshev approximation. Therefore, the residuals in the 
global field are evaluated for continuity equation <2.2. ?a> as 
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1 

(4.1 .2a) 

momentum 

equat ion 

(2.2.5b> as 
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1 < j < N-1 


(4.1 .2b) 


and at the boundaries 
upstream 



downstream 


(4.1 .2c> 


N 


= (X 


H 


- X ) h 

N-i N-2 


- (X - 
N 


X > 
N-2 


N-1 


+ (X 


N-1 


- X > 
N-2 


N 

(4.1 .2d> 


where h**** and are the guesses (from the finite difference 
preconditioner) of variables h and v for time level (p+1), h** and 
v’^ are the computed values at time p and d is a weighting factor 


which varies from 0.5 to 1. 



The root mean square <rms> value of the residual or 
function to be minimized is defined as the square root of 


the 


the 

inner product of the residuals obtained from (4.1.2) i.e. 



r H-1 


T 1/2 

r = 

1 + E 

fr . + r .1 

1 

rnns 

i N 

c j nr» 

J 


L j = l 

J 


The value of (4.1.5) gets minimized after some 
iterations and we get a more accurate solution to our governing 
equations. The whole process of iteration is discussed in the next 
sect ion. 


4* 2 Finit.e Difference Precondit-ioning 

As already mentioned, preconditioning is a good way to 
accelerate the convergence of iterative schemes designed to solve 
implicit spectral equations. The approximation matrix resulting 
from preconditioning should satisfy the following conditions. 

(a) the approximation matrix should be a fairly good 
representation of the original spectral operator 
matrix M 

ap 

(b) it should be easy and inexpensive to invert. 

The present work uses an implicit first-order finite 
difference preconditioner which employs forward-differences to 
discretize the spatial derivatives of the linearized governing 
equations. The linearization of the governing equation is carried 
out as follows. 

In the governing equation (2.2.3), the dependent 
variables h and v are replaced by (h + Ah) 


and 


( V + Av ) 
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respectively, where h and v are some approximation of h and v 
obtained from the finite difference implicit scheme and Ah and Av 
are the corrections that need to be made to h and v respectively. 
The equations are expanded in terms of corrections Ah and Av . The 
source term S is highly non-linear, and is therefore, linearized 

r 

by using Taylor's expansion. By dropping all second and higher 
order terms, the linearized equations become 


a Ah + b^Av h + b v Ah t c Ah v + c h Av 

ot OxOxOkOx 


R <4.2.1a> 
c 


d Av^ + e Av V + e v Av t f Ah - Ah (S ), - Av <S ) = - R 

OtOxOxOx rn rv M 


in which 


<4.2.1b> 


<S > 

r 


h 


ds 

"dh 



2 


V 


.7X3 

h 


(4. 2 .2a> 


and 


<S ) = 

r V 


dS 


-2 



V 

.4X3 

h 


<4.2 .2b) 


And and R^ are the residuals of the original equations for the 
approximate solution h and v. Subscripts x and t indicate the 
partial derivatives with respect to x and t respectively. 

If the residuals R and R are known for the current 

c M 

approximations h and v, we can use <4.1.2) to compute the 
corrections Ah and Av. We know, of course, only the spectral 
approximation of R^ and R^ which is computed as given above from 



h and v. However 
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, any attempt to solve (4.2.11 directly for the 
spectral approximation of Ah and Av will involve solutions of 
matrix equations that are full and ill-conditioned and 

computationally very expensive. The approach of finite difference 
preconditioning is to solve (4.2.1) by obtaining the values of Ah 
and Av at the collocation points by solving the finite difference 
version of (4.2.1) at each collocation points. 

Thus the right hand terms of equations (4.2.1) are 
computed by using (4.1.2) for the present time level values of h 
and V, and the left hand terms are discretized refer to finite 
difference form using Preissmann scheme again (But with a 

different weighting factor, fi) . After rearranging the 
corresponding terms, the final system of equations becooe 

a* Ah. + b° Av . +• c* Ah . + d* Av . = - R . (4.2.3a) 

J J J J i J + * i J + i C4 

1 < j < N-1 


a"* Ah + b” Av. + c”* Ah. + d’" Av. = - R . (4.2.3b) 

J j J J J J + A J J + 4 

1 < j < N-1 

and the equations at the boundaries are 


a'" Ah + b"" Av = - R (4.2.3c) 

1111 1 


a Ah 

N-2 N-2 


+ C 


N-1 


Ah 


N-1 


a^ Ah 

N N 


-R 


N 


(4.2. 3d) 


The marching techniques used to solve (4.2.3) generates 


at each marching step, a sparse matrix system of the form as 
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mentioned earlier C2-3.2.3). 

Again this matrix system is solved using a standard 
library subroutine and the corrections to the solution are thus 
evaluated. It must be noted that the use of insslicit finite 
difference schemes in the preconditioner ensures its unconditional 
stability. At the end of an iteration, we obtain the finite 
difference corrections at the collfxration points. Ah. and Av.. 

j i 

These corrections are multiplied by an acceleration 
parameter w and added to h and v to get the new guess of h and v. 
This entire process is repeated till the process converges, i.e. 
the rms of the residuals does not decrease any further. 

4.3 Preconditioned Residual Minimization 

The procedure adopted is as given in Chapter III, with 
the clarification that the residual R is made up of 

(1) the upstream boundary residual 

C2> the residuals of the continuity equation R . at the 

interior points 

(31 the residuals of the momentum equation R. . at the 

Mj 

interior points 

(4) the downstream boundary residual R . 

Therefore the inner product, say, (R, R ) is 

(R. R ) = R. r; - R„ R„ + V Rj r; 

J = i 


(4.3.1) 



3 



Fig. 4.1 


Flow Chart of the Alogrithm 



CHAPTER V 


RESULTS AND DISCUSSION 

In this Chapter, the four -point ioiplicit finite 
difference and the Chebyshev Collocation spectral methcKis are used 
to route a log pearson type III hydrograph through a wide 
rectangular open-channel. The nunerical results so obtained are 
assessed for their coRif»utational accuracy. Two hydrographs are 
used for this purpose. The peak discharge in both the hydrographs 
is equal to twenty tiokes the initial steady discharge. However, 
the peak discharge in the first hydrograph occurs at 72 hours 
while in the second it occurs at 20 hours. This makes the depth 
and velcHrity vary more rapidly in the second case as compared to 
in the first case. The time at the Centre of Cravity of the 
hydrograph is taken 1.2 tiBtes the time at the peak discharge 
tt =1.2) for both the cases. 

r 

The governing partial differential equations do not have 
an analytical solution for the above boundary condition and in the 
general case of friction affected flow. This makes it inupossible 
to exactly quantify the error introduced by the numerical scheme. 
Therefore, in the present study, the finite difference solution 
using 101 grid points and a ccmputational time step of lO minutes 
is taken as “exact solution". The numerical error in any other 
case is defined as below 


• 100 


% error = 


F - F* 

F^~' 


(5.1) 



36 


where^ F represents either the normalized depth or the normalized 
velocity and the subscript e indicates the "exact" result Equation 
(5. il is used to assess the resolution and tine stepping errors 
resulting fron the application of finite difference and spectral 
techniques for routing the flood wave. Three different values of 
the total number of grid points, N = 41, 21 and 11 are used for 
evaluating the spatial resolution errors- Two different values of 
the computational time step. At = 10 minutes and 4 hours are used 
for evaluating the time stepping errors. 

5.1 ^k]n-dinensionalizat ion of Equations: 

As mentioned, in Chapter II, St. Venant equations have 

been non -d i mens ionali zed to get the governing equations used in 

the present study. Using the parameters H , V and L representing 

o o o 

the initial uniform flow depth, initial uniform flow velocity and 
the channel length repectively, it is possible to obtain 12 
different sets of normalized equations. However, extensive 
numerical experimentation has shown that few of these sets give 
converging numerical results when used for routing the flood 
hydrograph. Further evaluation of numerical schemes is done only 
for the best set from the above which is selected by means of the 
following criteria : 

(a) Iteration procedure used for obtaining the spectral 
solution converges all times and for all ccMiputat ional 
time steps. 

(b) The rms value of residuals obtained from the finite 
difference solution is generally small. 
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Coefficients b,c,d,erf,9 and p in the 
ooooooo o 

normalized equations C2. X-'»> corresponding to all the 12 sets are 
given in Table 5.1. Remarks based on the earlier mentioned 
criteria for the "best set" are also included in this table. It 
can be seen that the set 4 is the "best set" and so it is chosen 
for further application. 

Table 5.1. Coefficients in Normalized St. Venant Equations 



1 

2 

3 

4 

5 

Set 

ao 

■Oi 

0 

U 

do 

eo 

H 

Ho 


Ho 

Vo 

vl 


VoTo 

IBH 

Lo 

g To 

g Lo 


Ho 

Ho 

Ho 

1 A 

VoTo 


VoTo 

Lo 

Lo 

JL a 

Lo 

■■H' 

Ho 

Ho 

Ho 

Vo Lo 

vl 


VoTo 

Lo 

Lo 

g Ho To 

g Ho 

M 

— 

VoTo 

VoTo 

Vo 

vl 


BBflii 

Lo 

Lo 

g To 

g Lo 

5 

o 

1 

VoTo 

Lo 

VoTo 

Lo 

mmm 

VoTo 

Lo 

A 

1 O 

VoTo 

VoTo 

r 

0 

< 

0 

vl 

o 

J. a W 

Lo 

Lo 

g Ho To 

g Ho 

T 

Lo 

1.0 


Vo 

vl 

# 

Vo To 

A a 

g To 

g Lo 


Lo 

1.0 


1 o 

VoTo 

d 

< 

0 

H 

0 

M. m 


Lo 

Q 

Lo 

o 

1 

tp-l 

1.0 

Lo Vo 

vl 

# 

Vo To 

g Ho To 

g Ho 

10 

Ho 

Ho 

Ho 

Lo 

1 o 

VoTo 

Lo 

Lo 

< 

0 

H 

0 

J> a W 

m 

1.0 

< 

0 

H 

0 

VoTo 

Lo 

1 o 

Lo 

Lo 

Vo To 

JL a V 

D 

|||_||||||_ 

1.0 



1.0 

1 Vo To 


Vo To 


(continued to next pagel 























































































38 


(continued from previous page> 
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5.2 Weighting Parameters a, /?, and 9 

The finite difference method presented in CJiapter II and 
the spectral method using the Chebysev Collocation technique 
presented in Chapters 111 and IV uses three vreighting parameters 
a, ft, and 9. These parameters are used in order to give a proper 
weightage to the solution towards the unknown time level (p+D- a 
is the weighting coefficient used while finding the initial 
solution by the Preissmann Scheme. 9 is the weighting coefficient 
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in the Chebysev Collocation technique used for finding the 
residuals of either the finite difference solution or the 
residuals during the iteration stage of spectral method. /? is the 
weighting coefficient in the finite difference preconditioning 
part of the spectral method. It should be noted that a = & = 0.5 
result in a second-order accuracy in time and a = 9 — 1.0 result 
in a ccMVPletely implicit first-order accurate method. Exhaustive 
numerical experimentation during the preliminary stage of BxiMfel 
development has indicated that a = 0.66, ft = 0.85 and 9 - 1.0 
result in gcsod numerical convergence and the lowest residuals. It 
proved to be very difficult to obtain convergent results with 9 = 
0.5. It can be seen from the above discussion that the finite 
difference method applied in the present study, is almost 
second-order accurate Cot — 0.66) while the spectral method is only 
first-order accurate i9 = 1.0) as far as time discretization is 
concerned. It should be remembered here that the spectral method 
as applied in the present study, achieves spectral accuracy only 
in space while it uses a finite difference time stepping. 

5.3 Results for equations without Source Term 

As mentioned in Qhapter 11, the right hand side of 
equation C2-2-3ttj) is termed as the Source term of the momentum 
equation. The governing partial difference equations are strictly 
hyperbolic only when this term is equal to zero- On the other 
hand, in many natural flow situations, the source term is more 
dominant than the other terms in the momentum equation. In such 
cases, it can be shown that the equations are of diffusion type 
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when the acceleration terms in the momentum equation are 
neglected. It can also be seen that the source term which includes 
the frictional effects tEg. 2. 2. 41 is highly non-linear. 
Therefore, difficulties are encountered while including the 
frictional effects in the numerical models. This is especially so, 
as experienced during the course of this study, idien one is 
interested in increasing the order of accuracy of numerical 
schemes. Therefore, the spectral technique for solving the St. 
Venant equations has been developed in two stages. In the first 
stage, it has been assumed that the channel is almost horizontal 
and frictionless, results for tidiich arc discussed in the present 
section. In the second stage of development, an attempt has been 
made to include the frictional and channel slope effects into the 
model. Results for the above case are discussed in the next 
section. 

The numerical results, i.e. the depth and velocity 
values, are obtained at all the grid points as functions of time. 
However, only the results for channel mid -sect ion are presented 
below. Similar trend is cAiserved at all the points. Figures 
5.1a and 5.1b show the normalized depth and normalized velocity as 
function of normalized time at the channel mid -sect ion for the 
first hydrograph. These results, obtained using the F.D scheme 
with 101 grid points and At = 10 minutes, are taken as "exact 
solutions”. In the same way. Figures 5.2a and 5.2b show the 
"exact solutions" for depth and velocity at the channel 
mid-section for the second hydrograph. It can be seen that 
although the peak depth is approximately the same in both cases. 



MOR. '/BLOCITY (v/Vo) > NOa DEPTH (h/Ho) 





NOR. VELOCITY (vA'O) > NOR. DEP^TH (h/Ho) 
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Fig. 5.2 "Exact Solution" for Hydrograph-2, At=10 minutes 
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it is attained at a much faster rate in the case of second 
hydrograph. 

In order to study the behaviour of spatial resolution 
errors, CDnw>uter runs were made using the same At (= 10 minutes) 
as before, but with a reduced number of grid points, 14 = 41, 21 
and 11. Figures 5.3a and 5.3b show the % error in depth 
computation tcalculated using Eg. 5.1) as a function of time for 
spectral and finite difference methods respectively. It is obvious 
from these figures that in both methods the results converge to 
"exact solution" as the nunriser of grid points are increased. 
However, looking the other way round, it can be observed that the 
finite^ difference solution CFig- 5.3b) deteriorates faster than 
the spectral solution CFig 5.3a) as the number of grid points is 
decreased. This is expected since spectral techniques 
theoretically give "infinite-order" accuracy and in the present 
case, we are minimizing the time stepping errors by taking a very 
low value of At. Similar trend is observed when 7 . error in 
velocity computation, for first hydrograph, is plotted (Figs. 5.4a 
and 5.4b) as a function of time. 

Figures 5.5a and 5.5b show the 7 . error in depth 
computation as a function of time for spectral and finite 
difference methods respectively for the second hydrograph. Here, 
again the finite difference solution, in general, deteriorates 
faster than the spectral solution as the number of nodes CN) is 
decreased. Hewever, just before the occurrence of peak depth, the 
error in the i^ectral solution t4.5X) is slightly more than the 
error in the finite difference solution (3.2X). This is due to the 
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(a) Spectral Solution 



(b) Finite Difference Solution 

Fig. 5.4 Spatial Resolution Error in Velocity Computation, 
Hydrograph-1 , At=10 minutes 
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(b) Finite Difference Solution 

Fig. 5.5 Spatial Resolution Error in Depth Computation, 
Hyclrograph-2 , At=10 minutes 
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fact that the finite difference scheme is almost second -order 
accurate in time while the spectral method is first-order 
accurate. It should be noted that the second hydrograph represents 
a faster rising flood and therefore, the time stepping errors are 
more predominant. 

The domenance of time-stepping errors over the spatial 
resolution errors becomes mare accentuated when the computational 
time step is increased. Figures 5.6a and 5.6b shew the error 
variation in the depth computation for the first hydrograph when 
the computational time step. At equals to four hours. It should be 
mentioned that this time step, though large is not unrealistically 
so, for the given hydrographs. The error in the spectral solution 


(Fig 5.6a) 

is in the range of 2% whereas the error 

in 

the 

finite 

difference 

solution (Fig 

5.6b) is only about 

1%. 

This 

trend 

becomes more clear «4ien 

computations are Btade 

for 

the 

faster 


rising second hydrograph. The error in the spectral solution is as 
much as 15% (Fig 5.7a) while the error in the finite difference 
solution (Fig 5.7b) is only in the range of 5%. This, clearly, 
indicates the dominance of time stepping errors for the cases 
studied here. 

5.4 Results for Equations with Source Term 

In this section, the results obtained from the 
application of ccsmplete St. Venant equations including the source 
terms for routing a flexsd hydrograph are presented. Results for 
only the first hydrograph are discussed. Similar trend was 
observed when computations were made for the second hydrograph. 
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Fig 5-S shows the normalized depth at the channel mid-section 
obtained using the Preissmann scheme with 101 grid points and At 
equal to 10 minutes. This is taken as the “exact solution" for 
further comparisons. It should be mentioned that the hydrograph 
used in this study is exactly the same as the one used by Ctxiley 
and Moin C19761 in their finite -element model. The solution is 
shown in Fig 5.8 is essentially the same as cdttained by them. 

Figures 5.9a and 5.9b show the percentage error in depth 
conwputation for spectral and finite difference methods, 
respectively, when the nuHriser of grid points is reduced. Three 
different N values equal to 41, 21 and 11 are used for studying 
the spatial resolution errors. It can be cAiserved that error in 
peak depth estimation (normalized time = 4.51 is insignificant 
even when as few as 11 nodes are used. This is true for both 
spectral and finite difference solutions. The numerical error 
occurs only in the rising limb of the hydrograph- It is observed 
that the errors in the spectral solutions t^ig 5.9a) are of the 
saiMB order of magnitude as the errors in finite difference 
solutions (Fig 5-9b)- It should be noted here that the frictional 
effects inhibit the free movement of hydrograph down the channel, 
thereby resulting in higher flow depths. For example, the peak 
normalized depth is approximately equal to 2.3 when the source 
terms are not included and it is approximately equal to 6 when the 
source terms arc included. However, the tinm of occurrence of of 
peak depth is almost saom* in both the cases. This indicates that 
the time variation in depth is more pronounced when the samB 
source terms arc included. As a consequence, the time-stepping 
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procedure becomes more in«>ortant than the spatial resolutions. 
Therefore, the spectral method which is only first-order accurate 
in time does not give significantly better results than those 
obtained using the second-order time-accurate finite difference 
scheme. This conclusion is reinforced by the results obtained for 
the second hydrograph (not shown here) . 

The dominance of time stepping errors is demonstrated 
more dramatically when the ccumputat ional time step. At = 4 hours. 
Figure 5.10 shows the error in depth computation as a function of 
time when At = 4 hours. As before, error in the peak depth 
estimation (normalized time = 4.5) is insignificant for both 
spectral and finite difference solutions. However, the error in 
the rising limb of the hydrograph is as much as 39)( when spectral 
method is used as compared to an error of 15% in the finite 
difference solution. The dcminance of time-stepping errors is 
clearly demonstrated by the fact that there is no ilmsirovement in 
the solution even when the number of grid points is increased from 
N = 11 to N = 41. 

5.5 Computational Time 

The total number of operations per time-step envolved in 
finite difference method and in spectral method is approximately 
equal to 30 N*+168 N and AO N*+350 N respectively. 

The ratio of the C.P.U. time per iteration taken by the 
finite difference method to that taken by the spectral iMithod is 


approximately 4.0 
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5.6 Conclusions 

From the above discussion, it can be concluded that a 
spectral method which achieves spectral accuracy only in Space is 
not good for solving unsteady shallow-water flow equations. A 
finite difference scheme which is only first-order accurate in 
space, but second-order accurate in time gives much better results 
with less amount of computational effort. Therefore, future work 
in this area should concentrate on developing fully spectral 
methods vrtiich achieve spectral accuracy in time also. 



CHAPTER VI 


SUMMARY 

This investigation presents a ^ectral Method for 
solving the one— dimensional shallow-water wave equations. The 
spectral method is based on the Oiebyshev Collocation technique 
and a first-order finite difference time stepping. The method 
developed in this study is applied to route a log-pearson type III 
hydrograph through a wide rectangular channel. The spectral 
solutions are compared with the solutions obtained using the 
Preissmann finite difference scheme which is first-order accurate 
in space and second-order accurate in time. An error analysis is 
0 iade through a parametric study. Following are the conclusion 
obtained from this study: 

(1) Both the spectral and the finite difference schemes do 
not introduce significant errors as far as peak depth 
simulation is concerned. Hcnurever, the numerical results 
show significant errors in the rising limb of the 
hydrograph. This is especially true tdien large 
computational time step and few grid points are used. 

t2> Finite difference solution deteriorates faster than the 
spectral solution vrtien the number of grid points is 
reduced while keeping the computational time step SHkall 
i.e. keeping time stepping errors low. 

<3> The spectral solution is worse than the finite 


difference solution whenever the time stepping errors 
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are predominant aswhen the hydrcsgraph has a fast rising 
limb. Time stepping errors are also significant when the 
effect of friction is predominant. 

It is suggested that the future work in this area should 
concentrate on developing fully spectral methods which achieve 
spectral accuracy in both space and time. The present model may 
also be modified so that it can be applied for diffusion routing 
in a more efficient manner. A significant application for the 
spectral methods would be in the area of sediment transport 
modeling where aggradation and degradation of channel beds occur 
over long reaches and over long periods of time. 
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APPENDIX A 


Details of equation (2.3.2. 1>; for convenience the bar 
has been drcw>ped from the barred variables without changing the 
nature of the variables. 


Equation (2.3.2. lal; 
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APPENDIX B 


Details of matrix equation (2- 3- 2. 3); 

for convenience the bar has been dropped from thi 
barred variables without changing the nature of the variables. 
Elements of the coefficient matrix; 
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